Variation in vocal behavior under varying stress levels
Marcelo Araya-Salas, PhD & Timothy Wright, PhD
2023-02-14
Source code and data found at https://github.com/maRce10/budgie_inbre_stress_vocal_output
Load packages
## add 'developer' to packages to be installed from github
x <- c("remotes", "lubridate", "readxl", "pbapply", "viridis", "ggplot2",
"kableExtra", "knitr", "formatR", "MASS", "sp", "GGally", "brms",
"lme4", "dplyr", "purrr", "forcats", "tidyr", "modelr", "tidybayes",
"cowplot", "ggrepel", "posterior", "ggridges", "maRce10/PhenotypeSpace")
source("~/Dropbox/R_package_testing/sketchy/R/load_packages.R")
load_packages(x)
source("~/Dropbox/R_package_testing/brmsish/R/html_summary.R")
source("~/Dropbox/R_package_testing/brmsish/R/helpers.R")
source("~/Dropbox/R_package_testing/brmsish/R/read_summary.R")Functions and global parameters
opts_knit$set(root.dir = "..")
# set evaluation false
opts_chunk$set(fig.width = 10, fig.height = 6, warning = FALSE, message = FALSE,
tidy = TRUE)
read_excel_df <- function(...) data.frame(read_excel(...))
# for reading months in english format
sl <- Sys.setlocale(locale = "en_US.UTF-8")
standard_error <- function(x) sd(x)/sqrt(length(x))
cols <- viridis(10, alpha = 0.7)
# print results brms models
summary_brm_model <- function(x, gsub.pattern = NULL, gsub.replacement = NULL,
xlab = "Effect size") {
summ <- summary(x)$fixed
fit <- x$fit
betas <- grep("^b_", names(fit@sim$samples[[1]]), value = TRUE)
hdis <- t(sapply(betas, function(y) hdi(fit@sim$samples[[1]][[y]])))
summ_table <- data.frame(summ, hdis, iterations = attr(fit, "stan_args")[[1]]$iter,
chains = length(attr(fit, "stan_args")))
summ_table <- summ_table[rownames(summ_table) != "Intercept",
c("Estimate", "Rhat", "Bulk_ESS", "l.95..CI", "u.95..CI",
"iterations", "chains")]
summ_table <- as.data.frame(summ_table)
summ_table$Rhat <- round(summ_table$Rhat, digits = 3)
summ_table$CI_low <- round(unlist(summ_table$l.95..CI), digits = 3)
summ_table$CI_high <- round(unlist(summ_table$u.95..CI), digits = 3)
summ_table$l.95..CI <- summ_table$u.95..CI <- NULL
out <- lapply(betas, function(y) data.frame(variable = y, value = sort(fit@sim$samples[[1]][[y]],
decreasing = FALSE)))
posteriors <- do.call(rbind, out)
posteriors <- posteriors[posteriors$variable != "b_Intercept",
]
if (!is.null(gsub.pattern) & !is.null(gsub.replacement))
posteriors$variable <- gsub(pattern = gsub.pattern, replacement = gsub.replacement,
posteriors$variable)
gg <- ggplot(data = posteriors, aes(y = variable, x = value, fill = stat(quantile))) +
geom_density_ridges_gradient(quantile_lines = TRUE, quantile_fun = HDInterval::hdi,
vline_linetype = 2) + theme_classic() + labs(x = xlab,
y = "Predictor") + xlim(range(c(min(summ_table[, "CI_low"]),
max(summ_table[, "CI_high"])), 0)) + geom_vline(xintercept = 0,
col = "red") + scale_fill_manual(values = c("transparent",
"lightblue", "transparent"), guide = "none") + ggtitle(x$formula)
if (!is.null(gsub.pattern) & !is.null(gsub.replacement))
rownames(summ_table) <- gsub(pattern = gsub.pattern, replacement = gsub.replacement,
rownames(summ_table))
summ_table$Rhat <- ifelse(summ_table$Rhat > 1.05, cell_spec(summ_table$Rhat,
"html", color = "white", background = "red", bold = TRUE,
font_size = 12), cell_spec(summ_table$Rhat, "html"))
signif <- summ_table[, "CI_low"] * summ_table[, "CI_high"] > 0
df1 <- kbl(summ_table, row.names = TRUE, escape = FALSE, format = "html",
digits = 3)
df1 <- row_spec(kable_input = df1, row = which(summ_table$CI_low *
summ_table$CI_high > 0), background = adjustcolor(cols[9],
alpha.f = 0.3))
df1 <- kable_styling(df1, bootstrap_options = c("striped", "hover",
"condensed", "responsive"), full_width = FALSE, font_size = 12)
print(df1)
print(gg)
}# read ext sel tab calls
sels <- read.csv("./data/processed/tailored_budgie_calls_sel_tab.csv")
# keep only spectrographic parameters
sels <- sels[, c("sound.files", "selec", "duration", "meanfreq", "sd",
"freq.median", "freq.IQR", "time.IQR", "skew", "kurt", "sp.ent",
"time.ent", "entropy", "meandom", "mindom", "maxdom", "dfrange",
"modindx", "meanpeakf")]
sels$ID <- sapply(sels$sound.files, function(x) strsplit(x, "_")[[1]][1])
sels$month <- sapply(sels$sound.files, function(x) strsplit(x, "_")[[1]][2])
sels$day <- sapply(sels$sound.files, function(x) strsplit(x, "_")[[1]][3])
sels$year <- sapply(sels$sound.files, function(x) strsplit(x, "_")[[1]][4])
sels$date <- paste(sels$day, substr(sels$month, 0, 3), sels$year,
sep = "-")
sels$date <- as.Date(sels$date, format = "%d-%b-%Y")
# acoustic measurements
areas_by_week <- readRDS("./data/processed/acoustic_space_area_by_individual_and_week.RDS")
indiv_ovlp <- readRDS("./data/processed/acoustic_space_density_overlap_to_first_week_by_individual.RDS")
indiv_ovlp$treatment <- factor(indiv_ovlp$treatment, levels = c("Control",
"Medium Stress", "High Stress"))
group_ovlp <- readRDS("./data/processed/acoustic_space_density_overlap_to_group_by_week.RDS")
# group_ovlp <- do.call(rbind, lapply(unique(group_ovlp$ID),
# function(x) { X <- group_ovlp[group_ovlp$ID == x, ]
# X$overlap.to.group <- X$overlap.to.group -
# X$overlap.to.group[X$week ==
# min(as.numeric(as.character(X$week)))] X$distance.to.group <-
# X$distance.to.group - X$distance.to.group[X$week ==
# min(as.numeric(as.character(X$week)))] return(X) }))
group_ovlp$treatment <- factor(group_ovlp$treatment, levels = c("Control",
"Medium Stress", "High Stress"))Counts per individual
df <- as.data.frame(table(sels$ID))
names(df) <- c("ID", "Sample_size")
df <- df[order(df$Sample_size, decreasing = FALSE), ]
kb <- kable(df, row.names = FALSE)
kb <- kable_styling(kb, bootstrap_options = c("striped", "hover",
"condensed", "responsive"))
print(kb)| ID | Sample_size |
|---|---|
| 132YGMM | 6 |
| 125YGMM | 12 |
| 160YGHM | 16 |
| 310YGHM | 21 |
| 393YGHM | 25 |
| 303YGHM | 37 |
| 398WBHM | 41 |
| 365YGHM | 44 |
| 399YGLM | 46 |
| 300YGHM | 47 |
| 270YGHM | 64 |
| 407YGHM | 97 |
| 386WBHM | 100 |
| 377WWLM | 108 |
| 367WWNM | 119 |
| 371YYLM | 149 |
| 397YGHM | 175 |
| 378YYLM | 195 |
| 404WBHM | 196 |
| 258YGHM | 223 |
| 389WWLM | 230 |
| 262YGHM | 306 |
| 400YGHM | 360 |
| 388YGLM | 377 |
| 327YYLM | 404 |
| 396YBHM | 455 |
| 403WBHM | 566 |
| 356WBHM | 761 |
| 361YGHM | 767 |
| 402YGHM | 849 |
| 395WBHM | 854 |
| 360YGHM | 900 |
| 390WBHM | 935 |
| 405YBHM | 975 |
| 385YBHM | 1256 |
| 394YBHM | 1693 |
Add metadata
metadat <- read_excel_df("./data/raw/INBREStress_MasterDataSheet_14Nov19.xlsx")
# head(metadat)
sels$ID[sels$ID == "125YGMM"] <- "125YGHM"
sels$ID[sels$ID == "394YBHM"] <- "394WBHM"
# setdiff(sels$ID, metadat$Bird.ID) setdiff(metadat$Bird.ID,
# sels$ID)
sels$treatment <- sapply(1:nrow(sels), function(x) {
metadat$Treatment[metadat$Bird.ID == sels$ID[x]][1]
})
sels$treatment.room <- sapply(1:nrow(sels), function(x) {
metadat$Treatment.Room[metadat$Bird.ID == sels$ID[x]][1]
})
sels$round <- sapply(1:nrow(sels), function(x) {
metadat$Round[metadat$Bird.ID == sels$ID[x]][1]
})
sels$source.room <- sapply(1:nrow(sels), function(x) {
metadat$Source.Room[metadat$Bird.ID == sels$ID[x]][1]
})
sels$record.group <- sapply(1:nrow(sels), function(x) {
metadat$Record.Group[metadat$Bird.ID == sels$ID[x]][1]
})
# add week
out <- lapply(unique(sels$round), function(x) {
Y <- sels[sels$round == x, ]
min_date <- min(Y$date)
week_limits <- min_date + seq(0, 100, by = 7)
Y$week <- NA
for (i in 2:length(week_limits)) Y$week[Y$date >= week_limits[i -
1] & Y$date < week_limits[i]] <- i - 1
return(Y)
})
sels <- do.call(rbind, out)
sels$cort.baseline <- sapply(1:nrow(sels), function(x) {
if (sels$week[x] == 1)
out <- metadat$D3.CORT.Baseline[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 2)
out <- metadat$D7.CORT.Baseline[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 3)
out <- metadat$D14.CORT.Baseline[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 4)
out <- metadat$D21.CORT.Baseline[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 5)
out <- metadat$D28.CORT.Baseline[metadat$Bird.ID == sels$ID[x]][1]
return(out)
})
sels$cort.stress <- sapply(1:nrow(sels), function(x) {
if (sels$week[x] == 1)
out <- metadat$D3.CORT.Stress[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 2)
out <- metadat$D7.CORT.Stress[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 3)
out <- metadat$D14.CORT.Stress[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 4)
out <- metadat$D21.CORT.Stress[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 5)
out <- metadat$D28.CORT.Stress[metadat$Bird.ID == sels$ID[x]][1]
return(out)
})
sels$stress.response <- sels$cort.stress - sels$cort.baseline
sels$weight <- sapply(1:nrow(sels), function(x) {
if (sels$week[x] == 1)
out <- metadat$D3.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 2)
out <- metadat$D7.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 3)
out <- metadat$D14.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 4)
out <- metadat$D21.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 5)
out <- metadat$D28.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]
return(out)
})
sels$breath.count <- sapply(1:nrow(sels), function(x) {
if (sels$week[x] == 1)
out <- metadat$D3.Breath.Count[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 2)
out <- metadat$D7.Breath.Count[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 3)
out <- metadat$D14.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 4)
out <- metadat$D21.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 5)
out <- metadat$D28.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]
return(out)
})
# remove week 5
sels <- sels[sels$week != 5, ]agg_dat <- aggregate(selec ~ ID + week, data = sels, length)
# compare to week 1 agg_dat$call.count <-
# sapply(1:nrow(agg_dat), function(x) { baseline <-
# agg_dat$selec[agg_dat$week == 1 & agg_dat$ID == agg_dat$ID[x]]
# if (length(baseline) > 0) change <- agg_dat$selec[x] -
# baseline else change <- agg_dat$selec[x] return(change) } )
# without comparing to week 1
agg_dat$call.count <- sapply(1:nrow(agg_dat), function(x) agg_dat$selec[x])
agg_dat$selec <- NULL
agg_dat$baseline.CORT <- sapply(1:nrow(agg_dat), function(x) {
baseline <- sels$cort.baseline[sels$week == 1 & sels$ID == agg_dat$ID[x]]
current <- sels$cort.baseline[sels$week == agg_dat$week[x] & sels$ID ==
agg_dat$ID[x]]
if (length(baseline) > 0 & length(current) > 0)
change <- mean(current) - mean(baseline) else change <- NA
return(change)
})
agg_dat$stress.response <- sapply(1:nrow(agg_dat), function(x) {
# baseline <- sels$stress.response[sels$week == 1 & sels$ID
# == agg_dat$ID[x]]
current <- sels$stress.response[sels$week == agg_dat$week[x] &
sels$ID == agg_dat$ID[x]]
# if (length(baseline) > 0 & length(current) > 0) change <-
# mean(current) - mean(baseline) else change <- NA
return(mean(current))
})
agg_dat$stress.CORT <- sapply(1:nrow(agg_dat), function(x) {
baseline <- sels$cort.stress[sels$week == 1 & sels$ID == agg_dat$ID[x]]
current <- sels$cort.stress[sels$week == agg_dat$week[x] & sels$ID ==
agg_dat$ID[x]]
if (length(baseline) > 0 & length(current) > 0)
change <- mean(current) - mean(baseline) else change <- NA
return(change)
})
agg_dat$weight <- sapply(1:nrow(agg_dat), function(x) {
baseline <- sels$weight[sels$week == 1 & sels$ID == agg_dat$ID[x]]
current <- sels$weight[sels$week == agg_dat$week[x] & sels$ID ==
agg_dat$ID[x]]
if (length(baseline) > 0 & length(current) > 0)
change <- mean(current) - mean(baseline) else change <- NA
return(change)
})
agg_dat$breath.rate <- sapply(1:nrow(agg_dat), function(x) {
baseline <- sels$breath.count[sels$week == 1 & sels$ID == agg_dat$ID[x]]
current <- sels$breath.count[sels$week == agg_dat$week[x] & sels$ID ==
agg_dat$ID[x]]
if (length(baseline) > 0 & length(current) > 0)
change <- mean(current) - mean(baseline) else change <- NA
return(change)
})
agg_dat$acoustic.area <- sapply(1:nrow(agg_dat), function(x) {
area <- areas_by_week$raref.area[areas_by_week$ID == agg_dat$ID[x] &
areas_by_week$week == agg_dat$week[x]]
if (length(area) < 1)
area <- NA
return(area)
})
agg_dat$acoustic.distance <- sapply(1:nrow(agg_dat), function(x) {
distance <- indiv_ovlp$distance.to.first.week[indiv_ovlp$ID ==
agg_dat$ID[x] & indiv_ovlp$week == agg_dat$week[x]]
if (length(distance) < 1)
distance <- NA
return(distance)
})
agg_dat$acoustic.overlap <- sapply(1:nrow(agg_dat), function(x) {
overlap <- indiv_ovlp$overlap.to.first.week[indiv_ovlp$ID == agg_dat$ID[x] &
indiv_ovlp$week == agg_dat$week[x]]
if (length(overlap) < 1)
overlap <- NA
return(overlap)
})
agg_dat$acoustic.overlap.to.group <- sapply(1:nrow(agg_dat), function(x) {
overlap <- group_ovlp$overlap.to.group[group_ovlp$ID == agg_dat$ID[x] &
group_ovlp$week == agg_dat$week[x]]
if (length(overlap) < 1)
overlap <- NA
return(overlap)
})
agg_dat$treatment <- sapply(1:nrow(agg_dat), function(x) unique(sels$treatment[sels$ID ==
agg_dat$ID[x]]))
agg_dat$round <- sapply(1:nrow(agg_dat), function(x) unique(sels$round[sels$ID ==
agg_dat$ID[x]]))Physiological parameters
breath.count <- stack(metadat[, c("D3.Breath.Count", "D7.Breath.Count",
"D14.Breath.Count", "D21.Breath.Count", "D28.Breath.Count")])
weight <- stack(metadat[, c("D3.Bird.Weight..g.", "D7.Bird.Weight..g.",
"D14.Bird.Weight..g.", "D21.Bird.Weight..g.", "D28.Bird.Weight..g.")])
cort.stress <- stack(metadat[, c("D3.CORT.Stress", "D7.CORT.Stress",
"D14.CORT.Stress", "D21.CORT.Stress", "D28.CORT.Stress")])
cort.baseline <- stack(metadat[, c("D3.CORT.Baseline", "D7.CORT.Baseline",
"D14.CORT.Baseline", "D21.CORT.Baseline", "D28.CORT.Baseline")])
stress <- data.frame(metadat[, c("Bird.ID", "Treatment", "Round",
"Treatment.Room")], week = breath.count$ind, breath.count = breath.count$values,
weight = weight$values, cort.stress = cort.stress$values, cort.baseline = cort.baseline$values,
stress.response = cort.stress$values - cort.baseline$values)
# head(stress)
stress$week <- factor(sapply(strsplit(as.character(stress$week), "\\."),
"[[", 1), levels = c("D3", "D7", "D14", "D21", "D28"))
stress$day <- as.numeric(gsub("D", "", as.character(stress$week)))
stress$round <- paste("Round", stress$Round)
stress$treatment <- factor(stress$Treatment, levels = c("Control",
"Medium Stress", "High Stress"))
# remove 5th week
stress <- stress[stress$week != "D28", ]
stress_l <- lapply(stress$Bird.ID, function(x) {
X <- stress[stress$Bird.ID == x, ]
X$breath.count <- X$breath.count - X$breath.count[X$week == "D3"]
X$weight <- X$weight - X$weight[X$week == "D3"]
X$cort.stress <- X$cort.stress - X$cort.stress[X$week == "D3"]
X$cort.baseline <- X$cort.baseline - X$cort.baseline[X$week ==
"D3"]
X$stress.response <- X$stress.response #- X$stress.response[X$week == 'D3']
return(X)
})
stress <- do.call(rbind, stress_l)
agg_stress <- aggregate(cbind(breath.count, weight, stress.response,
cort.baseline) ~ week + treatment, stress, mean)
agg_stress_se <- aggregate(cbind(breath.count, weight, stress.response,
cort.baseline) ~ week + treatment, stress, standard_error)
names(agg_stress_se) <- paste(names(agg_stress_se), ".se", sep = "")
agg_stress <- cbind(agg_stress, agg_stress_se[, 3:6])
agg_stress$Week <- 1:4
bs <- 22
gg_breath.count <- ggplot(data = agg_stress, aes(x = Week, y = breath.count,
fill = treatment)) + geom_bar(stat = "identity") + geom_errorbar(aes(ymin = breath.count -
breath.count.se, ymax = breath.count + breath.count.se), width = 0.1) +
scale_fill_viridis_d(begin = 0.1, end = 0.9) + facet_wrap(~treatment,
ncol = 3, scale = "fixed") + labs(y = "Mean change in breath\nrate (breaths/min)",
x = "Week") + theme_classic(base_size = bs) + theme(legend.position = "none")
gg_weight <- ggplot(data = agg_stress, aes(x = Week, y = weight, fill = treatment)) +
geom_bar(stat = "identity") + geom_errorbar(aes(ymin = weight -
weight.se, ymax = weight + weight.se), width = 0.1) + scale_fill_viridis_d(begin = 0.1,
end = 0.9) + facet_wrap(~treatment, ncol = 3, scale = "fixed") +
labs(y = "Mean change in \nweight (grams)", x = "Week") + theme_classic(base_size = bs) +
theme(legend.position = "none")
gg_cort.baseline <- ggplot(data = agg_stress, aes(x = Week, y = cort.baseline,
fill = treatment)) + geom_bar(stat = "identity") + geom_errorbar(aes(ymin = cort.baseline -
cort.baseline.se, ymax = cort.baseline + cort.baseline.se), width = 0.1) +
scale_fill_viridis_d(begin = 0.1, end = 0.9) + facet_wrap(~treatment,
ncol = 3, scale = "fixed") + labs(y = "Mean change in\nbaseline CORT (ng/mL)",
x = "Week") + theme_classic(base_size = bs) + theme(legend.position = "none")
gg_stress.response <- ggplot(data = agg_stress, aes(x = Week, y = stress.response,
fill = treatment)) + geom_bar(stat = "identity") + geom_errorbar(aes(ymin = stress.response -
stress.response.se, ymax = stress.response + stress.response.se),
width = 0.1) + scale_fill_viridis_d(begin = 0.1, end = 0.9) +
facet_wrap(~treatment, ncol = 3, scale = "fixed") + labs(y = "Stress response \nCORT (ng/mL)",
x = "Week") + theme_classic(base_size = bs) + theme(legend.position = "none")
cowplot::plot_grid(gg_breath.count, gg_weight, gg_cort.baseline, gg_stress.response)
# cowplot::ggsave2(filename =
# './output/physiological_parameters_tru_time_70dpi.jpeg')
# cowplot::ggsave2(filename =
# './output/physiological_parameters_tru_time_300dpi.jpeg', dpi
# = 300)Stats
Models: Predicted physio measure ~ treatment + week (continuous) + IndRandom
Variables (Difference from Week 1): weight, BR, baseline CORT, Stress CORT, Stress Response
responses <- c("baseline.CORT", "stress.response", "stress.CORT",
"weight", "breath.rate")
predictors <- c("~ treatment + week + (1|ID) + (1|round)")
formulas <- expand.grid(responses = responses, predictors = predictors,
stringsAsFactors = FALSE)
vars_to_scale <- c(responses, "week")
# remove week 1
sub_agg_dat <- agg_dat[agg_dat$week != 1, ]
for (i in vars_to_scale) sub_agg_dat[, vars_to_scale] <- scale(sub_agg_dat[,
vars_to_scale])
physio_models <- lapply(1:nrow(formulas), function(x) {
sub_dat <- sub_agg_dat[!is.na(sub_agg_dat[names(sub_agg_dat) ==
formulas$responses[x]]), ]
sub_dat
mod <- try(brm(formula = paste(formulas$responses[x], formulas$predictors[x]),
iter = 20000, silent = 2, data = sub_dat, control = list(adapt_delta = 0.9),
chains = 4), silent = TRUE)
return(mod)
})
names(physio_models) <- paste(formulas$responses, formulas$predictors)
saveRDS(physio_models, "./data/processed/physiological_response_models.RDS")physio_models <- readRDS("./data/processed/physiological_response_models.RDS")
ggl <- list()
for (x in 1:length(physio_models)) {
ggl[[length(ggl) + 1]] <- html_summary(physio_models[[x]], gsub.pattern = "b_treatment|b_",
gsub.replacement = "", highlight = FALSE, remove.intercepts = TRUE,
model.name = names(physio_models)[x])
print(ggl[[length(ggl)]])
}baseline.CORT ~ treatment + week + (1|ID) + (1|round)
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | Intercept-student_t(3, 0.1, 2.5) sd-student_t(3, 0, 2.5) sigma-student_t(3, 0, 2.5) | baseline.CORT ~ treatment + week + (1 | ID) + (1 | round) | 20000 | 4 | 1 | 10000 | 852 | 0 | 4603.542 | 2994.716 | 905461639 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| HighStress | 0.885 | 0.125 | 1.656 | 1.001 | 5810.772 | 11998.851 |
| MediumStress | 0.195 | -0.632 | 1.009 | 1.001 | 4603.542 | 2994.716 |
| week | -0.012 | -0.137 | 0.111 | 1 | 18860.740 | 19591.347 |
stress.response ~ treatment + week + (1|ID) + (1|round)
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | Intercept-student_t(3, -0.1, 2.5) sd-student_t(3, 0, 2.5) sigma-student_t(3, 0, 2.5) | stress.response ~ treatment + week + (1 | ID) + (1 | round) | 20000 | 4 | 1 | 10000 | 291 | 0 | 17259.93 | 7182.645 | 1959704000 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| HighStress | -0.155 | -0.742 | 0.443 | 1 | 21183.71 | 24652.050 |
| MediumStress | -0.078 | -0.674 | 0.537 | 1 | 17259.93 | 7182.645 |
| week | -0.108 | -0.290 | 0.075 | 1 | 42051.84 | 25878.502 |
stress.CORT ~ treatment + week + (1|ID) + (1|round)
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | Intercept-student_t(3, 0.1, 2.5) sd-student_t(3, 0, 2.5) sigma-student_t(3, 0, 2.5) | stress.CORT ~ treatment + week + (1 | ID) + (1 | round) | 20000 | 4 | 1 | 10000 | 1382 | 0 | 1215.206 | 295.917 | 1998039246 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| HighStress | -0.252 | -1.032 | 0.534 | 1 | 6526.781 | 14743.384 |
| MediumStress | 0.095 | -0.790 | 0.949 | 1.003 | 4011.007 | 3739.939 |
| week | -0.080 | -0.195 | 0.031 | 1.003 | 1215.206 | 295.917 |
weight ~ treatment + week + (1|ID) + (1|round)
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | Intercept-student_t(3, 0.1, 2.5) sd-student_t(3, 0, 2.5) sigma-student_t(3, 0, 2.5) | weight ~ treatment + week + (1 | ID) + (1 | round) | 20000 | 4 | 1 | 10000 | 660 | 0 | 7094.293 | 12124.98 | 815174156 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| HighStress | -0.413 | -1.232 | 0.405 | 1.001 | 7094.293 | 13001.57 |
| MediumStress | -0.190 | -1.090 | 0.710 | 1.001 | 7498.298 | 12124.98 |
| week | -0.086 | -0.211 | 0.041 | 1.001 | 11368.583 | 23277.30 |
breath.rate ~ treatment + week + (1|ID) + (1|round)
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | Intercept-student_t(3, -0.4, 2.5) sd-student_t(3, 0, 2.5) sigma-student_t(3, 0, 2.5) | breath.rate ~ treatment + week + (1 | ID) + (1 | round) | 20000 | 4 | 1 | 10000 | 536 | 0 | 17035.19 | 9553.508 | 1346894033 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| HighStress | 0.219 | -0.171 | 0.617 | 1 | 17035.19 | 11812.759 |
| MediumStress | 0.110 | -0.317 | 0.544 | 1 | 18478.62 | 9553.508 |
| week | -0.771 | -0.916 | -0.628 | 1 | 35543.50 | 23097.254 |
names(ggl) <- names(physio_models)
Barplot and effect sizes side-by-side
cowplot::plot_grid(gg_breath.count + theme_classic(base_size = 10) +
theme(legend.position = "none"), ggl[[grep("breath", names(ggl))]])
cowplot::plot_grid(gg_weight + theme_classic(base_size = 10) + theme(legend.position = "none"),
ggl[[grep("weight", names(ggl))]])
cowplot::plot_grid(gg_cort.baseline + theme_classic(base_size = 10) +
theme(legend.position = "none"), ggl[[grep("base", names(ggl))]])
cowplot::plot_grid(gg_stress.response + theme_classic(base_size = 10) +
theme(legend.position = "none"), ggl[[grep("response", names(ggl))]])
Takeaways
Breath rate decreases gradually with time across after the first week
Stress response is higher in “high stress” birds compared to first week
Acoustic space projection
t-SNE
scale_param <- scale(sels[, c("duration", "meanfreq", "sd", "freq.median",
"freq.IQR", "time.IQR", "skew", "kurt", "sp.ent", "time.ent",
"entropy", "meandom", "mindom", "maxdom", "dfrange", "modindx",
"meanpeakf")])
tsne <- Rtsne(scale_param, dims = 2, perplexity = 30, verbose = FALSE,
max_iter = 5000)
saveRDS(tsne, "./data/processed/tsne_on_acoustic_parameters_jun_2021.RDS")tsne <- readRDS("./data/processed/tsne_on_acoustic_parameters_jun_2021.RDS")
Y <- as.data.frame(tsne$Y)
names(Y) <- c("TSNE1", "TSNE2")
sels <- data.frame(sels, Y)
sels$treatment <- factor(sels$treatment, levels = c("Control", "Medium Stress",
"High Stress"))
ggplot(sels, aes(x = TSNE1, y = TSNE2, col = as.factor(treatment))) +
geom_point() + labs(color = "Treatment") + scale_color_viridis_d(alpha = 0.4) +
theme_classic(base_size = 25) + guides(colour = guide_legend(override.aes = list(size = 10)))
Behavioral parameters
agg_call.count <- aggregate(cbind(call.count, acoustic.overlap.to.group) ~
week + treatment, agg_dat, mean)
agg_behav <- aggregate(cbind(acoustic.area, acoustic.overlap) ~ week +
treatment, agg_dat, mean)
agg_call.count_se <- aggregate(cbind(call.count, acoustic.overlap.to.group) ~
week + treatment, agg_dat, standard_error)
agg_behav_se <- aggregate(cbind(acoustic.area, acoustic.overlap) ~
week + treatment, agg_dat, standard_error)
agg_behav_se <- merge(agg_call.count_se, agg_behav_se, all = TRUE)
names(agg_behav_se) <- paste(names(agg_behav_se), ".se", sep = "")
agg_behav <- merge(agg_call.count, agg_behav, all = TRUE)
agg_behav <- cbind(agg_behav, agg_behav_se[, 3:6])
bs <- 22
agg_behav$treatment <- factor(agg_behav$treatment, levels = c("Control",
"Medium Stress", "High Stress"))
gg_call.count <- ggplot(data = agg_behav, aes(x = week, y = call.count,
fill = treatment)) + geom_bar(stat = "identity") + geom_errorbar(aes(ymin = call.count -
call.count.se, ymax = call.count + call.count.se), width = 0.1) +
scale_fill_viridis_d(begin = 0.1, end = 0.9) + facet_wrap(~treatment,
ncol = 3, scale = "fixed") + labs(y = "Vocal output", x = "Week") +
theme_classic(base_size = bs) + theme(legend.position = "none")
gg_acoustic.area <- ggplot(data = agg_behav, aes(x = week, y = acoustic.area,
fill = treatment)) + geom_bar(stat = "identity") + geom_errorbar(aes(ymin = acoustic.area -
acoustic.area.se, ymax = acoustic.area + acoustic.area.se), width = 0.1) +
scale_fill_viridis_d(begin = 0.1, end = 0.9) + facet_wrap(~treatment,
ncol = 3, scale = "fixed") + labs(y = "Change in vocal diversity",
x = "Week") + theme_classic(base_size = bs) + theme(legend.position = "none")
gg_acoustic.overlap <- ggplot(data = agg_behav, aes(x = week, y = acoustic.overlap,
fill = treatment)) + geom_bar(stat = "identity") + geom_errorbar(aes(ymin = acoustic.overlap -
acoustic.overlap.se, ymax = acoustic.overlap + acoustic.overlap.se),
width = 0.1) + scale_fill_viridis_d(begin = 0.1, end = 0.9) +
facet_wrap(~treatment, ncol = 3, scale = "fixed") + labs(y = "Vocal Plasticity",
x = "Week") + theme_classic(base_size = bs) + theme(legend.position = "none")
gg_acoustic.overlap.group <- ggplot(data = agg_behav, aes(x = week,
y = acoustic.overlap.to.group, fill = treatment)) + geom_bar(stat = "identity") +
geom_errorbar(aes(ymin = acoustic.overlap.to.group - acoustic.overlap.to.group.se,
ymax = acoustic.overlap.to.group + acoustic.overlap.to.group.se),
width = 0.1) + scale_fill_viridis_d(begin = 0.1, end = 0.9) +
facet_wrap(~treatment, ncol = 3, scale = "fixed") + labs(y = "Vocal convergence",
x = "Week") + theme_classic(base_size = bs) + theme(legend.position = "none")
cowplot::plot_grid(gg_call.count, gg_acoustic.area, gg_acoustic.overlap,
gg_acoustic.overlap.group)
# cowplot::ggsave2(filename =
# './output/vocal_parameters_tru_time_70dpi.jpeg')
# cowplot::ggsave2(filename =
# './output/vocal_parameters_tru_time_300dpi.jpeg', dpi = 300)Stats
Model: Predicted behavior ~ treatment + week (continuous) + IndRandom
Variables: # calls, Distance moved from self in first week, Overlap to original acoustic space, Match to group repertoire, Maybe overall size of acoustic space
Do as comparison to week one using rarefacted calls and kernel density
responses <- c("call.count", "acoustic.area", "acoustic.overlap",
"acoustic.overlap.to.group")
predictors <- c("~ treatment + week + (1|ID) + (1|round)")
formulas <- expand.grid(responses = responses, predictors = predictors,
stringsAsFactors = FALSE)
vars_to_scale <- c(responses, "week")
for (i in vars_to_scale) agg_dat[, vars_to_scale] <- scale(agg_dat[,
vars_to_scale])
behav_models <- lapply(1:nrow(formulas), function(x) {
sub_dat <- agg_dat[!is.na(agg_dat[names(agg_dat) == formulas$responses[x]]),
]
# remove week 1
if (!grepl("count|group", formulas$responses[x]))
sub_dat <- sub_dat[sub_dat$week != 1, ]
mod <- brm(formula = paste(formulas$responses[x], formulas$predictors[x]),
iter = 50000, silent = 2, data = sub_dat, control = list(adapt_delta = 0.9,
max_treedepth = 15), chains = 4, prior = c(prior(normal(0,
5), "b"), prior(normal(0, 10), "Intercept"), prior(student_t(3,
0, 10), "sd"), prior(student_t(3, 0, 10), "sigma")))
return(mod)
})
names(behav_models) <- paste(formulas$responses, formulas$predictors)
saveRDS(behav_models, "./data/processed/behavioral_response_models.RDS")behav_models <- readRDS("./data/processed/behavioral_response_models.RDS")
ggl <- list()
for (x in 1:length(behav_models)) {
ggl[[length(ggl) + 1]] <- html_summary(behav_models[[x]], gsub.pattern = "b_treatment|b_",
gsub.replacement = "", highlight = FALSE, remove.intercepts = TRUE,
model.name = names(behav_models)[x])
print(ggl[[length(ggl)]])
}call.count ~ treatment + week + (1|ID) + (1|round)
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 5) Intercept-normal(0, 20) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20) | call.count ~ treatment + week + (1 | ID) + (1 | round) | 50000 | 4 | 1 | 25000 | 1540 | 0 | 25959.37 | 35763.57 | 1897146692 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| HighStress | 0.787 | 0.217 | 1.354 | 1 | 28006.98 | 44077.55 |
| MediumStress | 0.396 | -0.183 | 0.967 | 1 | 25959.37 | 35763.57 |
| week | -0.058 | -0.202 | 0.085 | 1 | 66399.00 | 56720.49 |
acoustic.area ~ treatment + week + (1|ID) + (1|round)
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 5) Intercept-normal(0, 20) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20) | acoustic.area ~ treatment + week + (1 | ID) + (1 | round) | 50000 | 4 | 1 | 25000 | 3780 | 0 | 6759.291 | 5875.617 | 76261279 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| HighStress | -0.218 | -1.058 | 0.633 | 1.001 | 6759.291 | 5875.617 |
| MediumStress | -0.707 | -1.558 | 0.151 | 1 | 9861.602 | 20907.370 |
| week | 0.035 | -0.083 | 0.151 | 1 | 7041.084 | 9145.749 |
acoustic.overlap ~ treatment + week + (1|ID) + (1|round)
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 5) Intercept-normal(0, 20) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20) | acoustic.overlap ~ treatment + week + (1 | ID) + (1 | round) | 50000 | 4 | 1 | 25000 | 2656 | 0 | 21041.63 | 31360.35 | 2038732099 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| HighStress | 0.900 | 0.050 | 1.741 | 1 | 22453.07 | 37240.04 |
| MediumStress | 0.236 | -0.692 | 1.155 | 1 | 21041.63 | 31360.35 |
| week | -0.169 | -0.307 | -0.031 | 1 | 57619.53 | 58721.25 |
acoustic.overlap.to.group ~ treatment + week + (1|ID) + (1|round)
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 5) Intercept-normal(0, 20) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20) | acoustic.overlap.to.group ~ treatment + week + (1 | ID) + (1 | round) | 50000 | 4 | 1 | 25000 | 4289 | 0 | 707.655 | 171.241 | 450263308 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| HighStress | 0.345 | -0.303 | 1.016 | 1.001 | 11286.051 | 23797.575 |
| MediumStress | 0.361 | -0.290 | 1.038 | 1.001 | 21649.351 | 39958.132 |
| week | -0.220 | -0.386 | -0.032 | 1.006 | 707.655 | 171.241 |
names(ggl) <- names(behav_models)
Barplot and effect sizes side-by-side
col_pointrange <- viridis(10)[2]
gg_coeffs <- lapply(behav_models, function(x) {
vars <- grep("b_", posterior::variables(x), value = TRUE)
draws <- posterior::as_draws_array(x, variable = vars)
coef_table <- draw_summary(draws, variables = vars, probs = c(0.025,
0.975), robust = TRUE)
coef_table$predictor <- rownames(coef_table)
coef_table$predictor <- gsub("b_treatment|b_", "", coef_table$predictor)
coef_table$predictor <- gsub("Stress", " stress", coef_table$predictor)
coef_table$predictor <- gsub("week", "Week", coef_table$predictor)
coef_table <- coef_table[coef_table$predictor != "Intercept",
]
gg_coef <- ggplot2::ggplot(data = coef_table, aes(x = Estimate,
y = predictor)) + geom_vline(xintercept = 0, lty = 2) + ggplot2::geom_point(size = 4,
col = col_pointrange) + ggplot2::geom_errorbar(ggplot2::aes(xmin = `l-95% CI`,
xmax = `u-95% CI`), width = 0, col = col_pointrange) + ggplot2::theme_classic(base_size = bs) +
ggplot2::theme(axis.ticks.length = ggplot2::unit(0, "pt"),
plot.margin = ggplot2::margin(0, 0, 0, 0, "pt"), legend.position = "none",
strip.background = ggplot2::element_blank(), strip.text = ggplot2::element_blank()) +
ggplot2::labs(x = "Effect size", y = "")
return(gg_coef)
})
# cowplot::plot_grid(gg_call.count + theme_classic(base_size =
# 20) + theme(legend.position = 'none'),
# gg_coeffs[[grep('count', names(gg_coeffs))]], nrow = 2)
# cowplot::plot_grid(gg_acoustic.area + theme_classic(base_size
# = 20) + theme(legend.position = 'none'),
# gg_coeffs[[grep('area', names(gg_coeffs))]])
# cowplot::plot_grid(gg_acoustic.overlap +
# theme_classic(base_size = 20) + theme(legend.position =
# 'none'), gg_coeffs[[grep('overlap ~', names(gg_coeffs))]])
# cowplot::plot_grid(gg_acoustic.overlap.group +
# theme_classic(base_size = 10) + theme(legend.position =
# 'none'), gg_coeffs[[grep('group', names(gg_coeffs))]])
bs <- 10
cowplot::plot_grid(gg_call.count + theme_classic(base_size = bs) +
theme(legend.position = "none"), gg_acoustic.area + theme_classic(base_size = bs) +
theme(legend.position = "none"), gg_acoustic.overlap + theme_classic(base_size = bs) +
theme(legend.position = "none"), gg_acoustic.overlap.group + theme_classic(base_size = bs) +
theme(legend.position = "none"), gg_coeffs[[grep("count", names(gg_coeffs))]] +
theme_classic(base_size = bs), gg_coeffs[[grep("area", names(gg_coeffs))]] +
theme_classic(base_size = bs), gg_coeffs[[grep("overlap ~", names(gg_coeffs))]] +
theme_classic(base_size = bs), gg_coeffs[[grep("group", names(gg_coeffs))]] +
theme_classic(base_size = bs), nrow = 2, rel_heights = c(1.8,
1))
# cowplot::ggsave2(filename =
# './output/bar_graphs_and_estimates_70dpi.jpeg', width = 25,
# height = 9) cowplot::ggsave2(filename =
# './output/bar_graphs_and_estimates_300dpi.jpeg', dpi = 300,
# width = 25, height = 9)Takeaways
Higher vocal output in “high stress” birds compared to control
Higher acoustic overlap to themselves in week 1 for “high stress” birds compared to control
Decrease in overlap to themselves in week 1 across time
Session information
## R version 4.1.0 (2021-05-18)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.2 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/atlas/libblas.so.3.10.3
## LAPACK: /usr/lib/x86_64-linux-gnu/atlas/liblapack.so.3.10.3
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=pt_BR.UTF-8
## [7] LC_PAPER=es_CR.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=es_CR.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] PhenotypeSpace_0.1.0 warbleR_1.1.26 NatureSounds_1.0.4
## [4] seewave_2.1.8 tuneR_1.3.3.1 ggridges_0.5.3
## [7] posterior_1.2.0 ggrepel_0.9.1 cowplot_1.1.1
## [10] tidybayes_3.0.2 modelr_0.1.8 tidyr_1.2.0
## [13] forcats_0.5.1 purrr_0.3.4 dplyr_1.0.8
## [16] lme4_1.1-28 Matrix_1.3-4 brms_2.16.3
## [19] Rcpp_1.0.8 GGally_2.1.2 sp_1.4-6
## [22] MASS_7.3-54 formatR_1.11 knitr_1.37
## [25] kableExtra_1.3.4 ggplot2_3.3.5 viridis_0.6.2
## [28] viridisLite_0.4.0 pbapply_1.5-0 readxl_1.3.1
## [31] lubridate_1.8.0 remotes_2.4.2
##
## loaded via a namespace (and not attached):
## [1] utf8_1.2.2 tidyselect_1.1.2 fftw_1.0-6.1
## [4] htmlwidgets_1.5.4 grid_4.1.0 munsell_0.5.0
## [7] codetools_0.2-18 DT_0.20 miniUI_0.1.1.1
## [10] withr_2.4.3 spatstat.random_2.1-0 Brobdingnag_1.2-7
## [13] colorspace_2.0-3 highr_0.9 uuid_1.1-0
## [16] rstudioapi_0.13 stats4_4.1.0 dtw_1.22-3
## [19] tensor_1.5 bayesplot_1.8.1 labeling_0.4.2
## [22] emmeans_1.8.1-1 rstan_2.21.3 polyclip_1.10-0
## [25] farver_2.1.0 bridgesampling_1.1-2 rprojroot_2.0.2
## [28] coda_0.19-4 vctrs_0.3.8 generics_0.1.2
## [31] TH.data_1.1-0 xfun_0.29 R6_2.5.1
## [34] markdown_1.1 gamm4_0.2-6 projpred_2.0.2
## [37] bitops_1.0-7 spatstat.utils_2.3-0 reshape_0.8.8
## [40] assertthat_0.2.1 promises_1.2.0.1 scales_1.1.1
## [43] multcomp_1.4-17 gtable_0.3.0 goftest_1.2-3
## [46] processx_3.5.2 xaringanExtra_0.7.0 sandwich_3.0-1
## [49] rlang_1.0.1 systemfonts_1.0.4 splines_4.1.0
## [52] spatstat.geom_2.3-2 broom_0.7.12 checkmate_2.0.0
## [55] inline_0.3.19 yaml_2.3.5 reshape2_1.4.4
## [58] abind_1.4-5 threejs_0.3.3 crosstalk_1.2.0
## [61] backports_1.4.1 httpuv_1.6.5 rsconnect_0.8.25
## [64] tensorA_0.36.2 tools_4.1.0 ellipsis_0.3.2
## [67] spatstat.core_2.4-0 raster_3.5-15 jquerylib_0.1.4
## [70] RColorBrewer_1.1-2 proxy_0.4-26 plyr_1.8.6
## [73] base64enc_0.1-3 RCurl_1.98-1.6 ps_1.6.0
## [76] prettyunits_1.1.1 rpart_4.1-15 deldir_1.0-6
## [79] zoo_1.8-9 cluster_2.1.2 magrittr_2.0.2
## [82] ggdist_3.1.0 colourpicker_1.1.1 mvtnorm_1.1-3
## [85] matrixStats_0.61.0 shinyjs_2.1.0 mime_0.12
## [88] evaluate_0.15 arrayhelpers_1.1-0 xtable_1.8-4
## [91] shinystan_2.5.0 gridExtra_2.3 rstantools_2.1.1
## [94] compiler_4.1.0 tibble_3.1.6 crayon_1.5.0
## [97] minqa_1.2.4 StanHeaders_2.21.0-7 htmltools_0.5.2
## [100] mgcv_1.8-36 later_1.3.0 RcppParallel_5.1.5
## [103] DBI_1.1.1 boot_1.3-28 permute_0.9-7
## [106] cli_3.2.0 parallel_4.1.0 igraph_1.2.11
## [109] pkgconfig_2.0.3 signal_0.7-7 terra_1.5-21
## [112] spatstat.sparse_2.1-0 xml2_1.3.3 svUnit_1.0.6
## [115] dygraphs_1.1.1.6 svglite_2.1.0 bslib_0.3.1
## [118] webshot_0.5.2 estimability_1.4.1 rvest_1.0.2
## [121] stringr_1.4.0 distributional_0.3.0 callr_3.7.0
## [124] digest_0.6.29 vegan_2.5-7 spatstat.data_2.1-2
## [127] rmarkdown_2.11 cellranger_1.1.0 shiny_1.7.1
## [130] gtools_3.9.2 rjson_0.2.21 nloptr_2.0.0
## [133] lifecycle_1.0.1 nlme_3.1-152 jsonlite_1.8.0
## [136] fansi_1.0.2 pillar_1.7.0 lattice_0.20-44
## [139] loo_2.4.1 fastmap_1.1.0 httr_1.4.2
## [142] pkgbuild_1.3.1 survival_3.2-11 glue_1.6.1
## [145] xts_0.12.1 shinythemes_1.2.0 stringi_1.7.6
## [148] sass_0.4.0